#### Importing data and building core analysis data.frame

library(tidyverse)
library(readxl)
library(magrittr)
library(WDI) 
library(broom)
library(fixest)
library(ggthemes)
library(interflex)


select <- dplyr::select
shift <- data.table::shift

# setwd("rowan_extreme_weather_replication")


## Workflow
# 1. Load in all country-level outcomes
# 2. Load in all country-level stratifiers
# 3. Load in all country-level temperature data
# 4. Re-shape temperature data

#### Load in all country-level outcomes ####

## Grantham Climate Laws of the World
# https://www.climate-laws.org/legislation_and_policies.csv

# Load in the raw data from Grantham
grantham <- read_csv("data/grantham_climate_laws.csv") %>% 
  select(iso3c= `Geography ISO`, events = Events) %>% 
  mutate(date_1 = str_sub(events, start = 1L, end = 10L), # Take first date
         year_1 = str_sub(date_1, start = -4L, end = -1L), # Take first year
         amended = str_sub(events, start = 12L, end = -1L), # Keep information after first date
         date_2 = str_extract(amended, "\\d+/\\d+/\\d+"), # Take the next date from the amended half
         year_2 = str_sub(date_2, start = -4L, end = -1L), # Take the second year
         second_half = str_sub(amended, start = 34L, end = -1L), # Keep information after the second date
         date_3 = str_extract(second_half, "\\d+/\\d+/\\d+"), # Take the third date
         year_3 = str_sub(date_3, start = -4L, end = -1L)) %>% # Take the third year
  select(iso3c, starts_with("year")) %>% 
  mutate_at(vars(year_1, year_2, year_3), ~as.numeric(.))

# Extract IDs to make a panel
grantham_iso <- unique(grantham$iso3c)
grantham_years <- seq(from = min(grantham$year_1, na.rm=T), 
                      to = max(grantham[,2:4], na.rm=T),
                      by = 1)

# Combine into a panel
climate_laws <- full_join(
  # Create an empty panel
  cbind.data.frame( 
    iso3c = rep(grantham_iso,  
                each = length(grantham_years)),
    year = rep(grantham_years,
               times = length(grantham_iso))),
  # Merge the Grantham data in
  grantham, by = "iso3c") %>% 
  arrange(iso3c, year, year_1) %>% 
  rowwise() %>% 
  mutate(policy_1 = ifelse(year >= year_1, 1, 0), # Create binary for policy adoption
         policy_2 = ifelse(year >= year_2, 1, 0), # Create binary for policy amendment
         policy_3 = ifelse(year >= year_3, 1, 0), # Create binary for policy amendment
         policy_count = sum(policy_1, policy_2, policy_3, na.rm = T)) %>% # Count operating policies
  ungroup() %>% 
  group_by(iso3c, year) %>% 
  summarise(climate_laws = sum(policy_count, na.rm=T)) %>% # Summarise by year
  # Scale outcome variable 
  ungroup() %>% 
  mutate(scaled_climate_laws = scale(.$climate_laws)) %>%
  as.data.frame()

# Remove intermediate data
rm(grantham, grantham_iso, grantham_years)

## Load in climate policy data from Baldwin et al. 2019
# https://doi.org/10.1016/j.worlddev.2019.03.012

climate_policies <- read_csv("data/baldwin_wd_2019_policies.csv") %>% 
  mutate(scaled_fit = scale(.$fit),
         scaled_quota = scale(.$quota),
         scaled_carbon_priced = scale(.$carbon_priced)) %>%
  arrange(iso3c, year)

## Load in renewable electricity data
# 'EG.ELC.RNEW.ZS' from World Bank World Development Indicators

wdi_renew <- read_csv("data/wdi_renew.csv") %>% 
  mutate(scaled_elec_renew = scale(elec_renew)) %>% 
  select(-country)

## Load in carbon pricing data from World Bank's Carbon Pricing Dashboard (cf. Dolphin et al. 2020)

cp_list <- read_excel("data/wb_carbon_pricing_dashboard.xlsx", 
                      sheet = 1,
                      skip = 2) %>% 
  filter(`Type of juridiction covered`=="National") %>% 
  select(name = 1)

cp_price <- read_excel("data/wb_carbon_pricing_dashboard.xlsx", 
                       sheet = 3,
                       skip = 2) %>% 
  pivot_longer(names_to = "label", values_to = "price", `Price_label_1_1990`:`Price_rate_2_2020`) %>% 
  mutate(label2 = str_sub(label, start = 1L, end = -6L)) %>% 
  filter(label2 == "Price_rate_1") %>% 
  mutate(year = str_sub(label, start = -4L)) %>% 
  select(name = 1, year, price) %>% 
  mutate(year = as.numeric(year),
         price = as.numeric(price))

cp_scope <- read_excel("data/wb_carbon_pricing_dashboard.xlsx", 
                       sheet = 5,
                       skip = 2) %>% 
  select(name = 1, scope = 9) %>% 
  mutate(scope = str_sub(scope, start = 1L, end = -2L),
         scope = as.numeric(scope)/100)

ecp <- full_join(cp_price, cp_scope, by = c("name")) %>% 
  filter(name %in% cp_list$name) %>% 
  mutate(ecp = price*scope) %>% 
  mutate(name = str_remove(name, " carbon tax"),
         name = str_remove(name, " ETS"),
         name = str_remove(name, " carbon price floor"),
         name = str_remove(name, " ERF Safeguard Mechanism"),
         name = str_remove(name, " federal OBPS"),
         name = str_remove(name, " federal fuel charge"),
         name = str_remove(name, " pilot ETS")) %>%
  mutate(ecp_zeros = ifelse(is.na(ecp), 0, ecp)) %>% 
  mutate(iso3c = countrycode::countrycode(name, 
                                           origin = "country.name", 
                                           destination = "iso3c")) %>% 
  group_by(iso3c, year) %>% 
  summarise(ecp_zeros = sum(ecp_zeros)) %>% 
  ungroup() %>% 
  # Scale DV %>% 
  mutate(scaled_ecp_zeros = scale(.$ecp_zeros)) %>% 
  select(iso3c, year, ecp_zeros, scaled_ecp_zeros)
rm(cp_list, cp_price, cp_scope)

## Load in COP delegation data from Kaya and Schofield
# https://doi.org/10.1093/fpa/orz031

cop <- read_csv("data/kaya_fpa_2020.csv") %>% 
  select(iso3c = CountryCode, year = Year, delegates = Delegates0) %>% 
  arrange(iso3c, year) %>% 
  # Scale outcome variable 
  mutate(scaled_delegates = scale(.$delegates)) %>%
  as.data.frame()

## Load in climate institutional membership data from Rowan
# https://doi.org/10.1093/isq/sqaa092

gcg <- read_csv("data/rowan_isq_2020.csv") %>%
  select(iso3c, year, gcg = sum_climate_memberships) %>% 
  filter(year != "Cross-section") %>% 
  mutate(year = as.numeric(year)) %>% 
  arrange(iso3c, year) %>% 
  mutate(scaled_gcg = scale(.$gcg)) %>%
  as.data.frame()

## Load in climate finance data from OECD
# http://www.oecd.org/dac/financing-sustainable-development/development-finance-topics/climate-change.htm

finance_donor <- read_csv("data/climate_finance_donor.csv") %>% 
  mutate_at(vars(cf_multilateral_provided, cf_principal_provided, cf_significant_provided),
            replace_na, 0) %>% 
  mutate_at(vars(cf_multilateral_provided, cf_principal_provided, cf_significant_provided),
            round, 0) %>% 
  rowwise() %>% 
  mutate(cf_total_provided = sum(cf_multilateral_provided + cf_principal_provided + cf_significant_provided)) %>% 
  ungroup() %>% 
  # Log-transform
  mutate_at(vars(cf_principal_provided, cf_total_provided),
            list(~ log(1 + .))) %>% 
  select(iso3c, year, cf_principal_provided, cf_total_provided) %>% 
  distinct(iso3c, year, .keep_all = TRUE) %>%
  # Scale outcome variable 
  mutate(scaled_cf_total_provided = scale(.$cf_total_provided),
         scaled_cf_principal_provided = scale(.$cf_principal_provided)) %>%
  as.data.frame()

finance_recipient <- read_csv("data/climate_finance_recipient.csv") %>% 
    mutate_at(vars(cf_multilateral_received, cf_principal_received, cf_significant_received),
            replace_na, 0) %>% 
  mutate_at(vars(cf_multilateral_received, cf_principal_received, cf_significant_received),
            round, 0) %>% 
  rowwise() %>% 
  mutate(cf_total_received = sum(cf_multilateral_received + cf_principal_received + cf_significant_received)) %>% 
  ungroup() %>% 
  # Log-transform
  mutate_at(vars(cf_principal_received, cf_total_received),
            list(~ log(1 + .))) %>% 
  select(iso3c, year, cf_principal_received, cf_total_received) %>% 
  distinct(iso3c, year, .keep_all = TRUE) %>%
  # Scale outcome variable 
  mutate(scaled_cf_total_received = scale(.$cf_total_received),
         scaled_cf_principal_received = scale(.$cf_principal_received)) %>%
  as.data.frame()


#### Load in all country-level covariates/stratifiers ####  

rich <- read_csv("data/rich_dummy2.csv") %>% mutate(gdp_pcmkt = log(gdp_pcmkt))
democracy <- read_csv("data/democracy_dummy2.csv") %>% select(-year)
# agricultural <- read_csv("data/agricultural_dummy.csv")
# small <- read_csv("data/small_dummy.csv")
# vulnerable <- read_csv("data/vulnerable_dummy.csv")
# annex1 <- read_csv("data/kaya_fpa_2020.csv") %>% 
  # select(iso3c = CountryCode, annex1 = Annex1Dummy) %>% 
  # distinct()


#### Load in temperature regressors ####

## Make list of files to load in
# Change working folder
setwd("data/ms_daily_temp_pop_weighted")
myFiles <- list.files(pattern ="*.csv")

## Load in files 
temp_raw <- myFiles %>%
  map_df(~data.table::fread(.)) %>% 
  separate(date, into = c("year", "month", "day")) %>% 
  mutate_at(vars(year:day), list(~ as.numeric(.))) %>% 
  mutate(season = ifelse(month==12|month==1|month==2, "DJF",
                         ifelse(month==3|month==4|month==5, "MAM",
                                ifelse(month==6|month==7|month==8, "JJA", "SON")))) %>% 
  as.data.frame() 
# object.size(temp_raw) %>% format(., units = "Mb")

# Clean up ISO codes
temp_raw$iso3c <- countrycode::countrycode(sourcevar = temp_raw$iso,
                                           origin = "iso3c",
                                           destination = "iso3c")
# AU1, CH1, CYN, DN1, FI1, FR1, GB1, KAS, KOS, NL1, NZ1, SAH, SDS, SOL, US1

temp_raw$iso3c <- as.character(temp_raw$iso)
temp_raw$iso3c[temp_raw$iso=="AU1"] <- "AUS" # Australia
temp_raw$iso3c[temp_raw$iso=="CH1"] <- "CHN" # China
temp_raw$iso3c[temp_raw$iso=="DN1"] <- "DNK" # Denmark
temp_raw$iso3c[temp_raw$iso=="FI1"] <- "FIN" # Finland
temp_raw$iso3c[temp_raw$iso=="FR1"] <- "FRA" # France
temp_raw$iso3c[temp_raw$iso=="GB1"] <- "GBR" # Great Britain
temp_raw$iso3c[temp_raw$iso=="NL1"] <- "NLD" # Netherlands
temp_raw$iso3c[temp_raw$iso=="NZ1"] <- "NZL" # New Zealand
temp_raw$iso3c[temp_raw$iso=="US1"] <- "USA" # United States
temp_raw$iso3c[temp_raw$iso=="CYN"] <- NA  
temp_raw$iso3c[temp_raw$iso=="KAS"] <- NA  
temp_raw$iso3c[temp_raw$iso=="KOS"] <- NA  
temp_raw$iso3c[temp_raw$iso=="SAH"] <- NA  
temp_raw$iso3c[temp_raw$iso=="SDS"] <- NA  
temp_raw$iso3c[temp_raw$iso=="SOL"] <- NA  
temp_raw <- filter(temp_raw, !is.na(iso3c)) 

# Change working folder
setwd("../..")

#### Tidy the data to make the right regressors ####

#' Mean annual temp
#' Mean hottest season's seasonal temp
#' Degree days in hottest season
#' --> Doing this all in one pipe

# Identify hottest season within a country
hottest_season <- temp_raw %>% 
  group_by(iso3c, season) %>% 
  summarise(mean_season_temp = mean(temperature)) %>% 
  ungroup() %>%
  group_by(iso3c) %>% 
  mutate(mean_hottest_season_temp = max(mean_season_temp)) %>% 
  filter(mean_season_temp == mean_hottest_season_temp) %>% 
  mutate(hottest_season = season) %>% 
  select(iso3c, hottest_season, mean_hottest_season_temp)
# 200 x 3


# Re-load temperature in new shape for spatial terms

temp_df <- full_join(temp_raw, hottest_season, by = "iso3c") %>% 
  arrange(iso3c, year, month, day) %>% 
  select(iso3c, year, month, day, season, hottest_season, temperature) %>% 
  # Get annual average temperature
  group_by(iso3c, year) %>% 
  mutate(annual_average_temp = mean(temperature, na.rm=T)) %>% 
  ungroup() %>% 
  # Select only the days in the hottest season
  filter(season == hottest_season) %>% 
  # Temperature years run December--November for seasonality
  # Offset December to match seasonal year to year
  mutate(season_year = ifelse(month == 12, year + 1, year)) %>% 
  # Drop year and rename season_year for merging
  select(-year) %>% 
  rename(year = season_year) %>% 
  filter(year<2019) %>% 
  # Group by country and season_year
  group_by(iso3c, year) %>% 
  # Create the annual average seasonal temp in the hottest season
  mutate(hottest_season_temp = mean(temperature)) %>% 
  ungroup() %>% 
  # Create dummies for daily temperature realizations
  mutate(temp_bin_10 = ifelse(temperature <= 10, 1, 0),
         temp_bin_15 = ifelse(temperature > 10 & temperature <= 15, 1, 0),
         temp_bin_20 = ifelse(temperature > 15 & temperature <= 20, 1, 0),
         temp_bin_25 = ifelse(temperature > 20 & temperature <= 25, 1, 0),
         temp_bin_30 = ifelse(temperature > 25 & temperature <= 30, 1, 0),
         temp_bin_35 = ifelse(temperature > 30 & temperature <= 35, 1, 0),
         temp_bin_40 = ifelse(temperature > 35, 1, 0)) %>% 
  # Aggregate the number of days per year in each realized temperature bin
  group_by(iso3c, year) %>%
  mutate(count_temp_bin_10 = sum(temp_bin_10),
         count_temp_bin_15 = sum(temp_bin_15),
         count_temp_bin_20 = sum(temp_bin_20),
         count_temp_bin_25 = sum(temp_bin_25),
         count_temp_bin_30 = sum(temp_bin_30),
         count_temp_bin_35 = sum(temp_bin_35),
         count_temp_bin_40 = sum(temp_bin_40)) %>% 
  ungroup() %>% 
  # Keep only the important variables, tidy
  select(iso3c, year,
         annual_average_temp, hottest_season_temp, starts_with("count")) %>% 
  # Collapse to the year level, drop days and months
  distinct(iso3c, year, .keep_all = T) %>% 
  # Arrange for lags
  arrange(iso3c, year) %>% 
  group_by(iso3c) %>% 
  # Set lags
  do(data.frame(., setNames(shift(.$annual_average_temp, 1:3),
                            paste0('annual_average_temp_lag', 1:3)))) %>% 
  do(data.frame(., setNames(shift(.$hottest_season_temp, 1:3),
                            paste0('hottest_season_temp_lag', 1:3)))) %>% 
  do(data.frame(., setNames(shift(.$count_temp_bin_10, 1:3), paste0('count_temp_bin_10_lag', 1:3)))) %>%
  do(data.frame(., setNames(shift(.$count_temp_bin_15, 1:3), paste0('count_temp_bin_15_lag', 1:3)))) %>%
  do(data.frame(., setNames(shift(.$count_temp_bin_20, 1:3), paste0('count_temp_bin_20_lag', 1:3)))) %>%
  do(data.frame(., setNames(shift(.$count_temp_bin_25, 1:3), paste0('count_temp_bin_25_lag', 1:3)))) %>%
  do(data.frame(., setNames(shift(.$count_temp_bin_30, 1:3), paste0('count_temp_bin_30_lag', 1:3)))) %>%
  do(data.frame(., setNames(shift(.$count_temp_bin_40, 1:3), paste0('count_temp_bin_35_lag', 1:3)))) %>%
  do(data.frame(., setNames(shift(.$count_temp_bin_35, 1:3), paste0('count_temp_bin_40_lag', 1:3)))) %>% 
  ungroup() %>% 
  as.data.frame() %>% 
  # Group by country  
  group_by(iso3c) %>%
  # Standardize temperature within countries
  mutate(z_annual_average_temp = scale(annual_average_temp),
         z_hottest_season_temp = scale(hottest_season_temp)) %>%
  # Ungroup, keep only observation period, and recode heat shock as positive only
  ungroup() %>% 
  filter(year >= 1990) %>% 
  mutate(annual_heat_shock = ifelse(z_annual_average_temp > 0, z_annual_average_temp, 0),
         seasonal_heat_shock = ifelse(z_hottest_season_temp > 0, z_hottest_season_temp, 0)) %>% 
  # Sum the heat shocks and create a lag
  group_by(iso3c) %>%
  mutate(annual_heat_shock_cumsum = cumsum(annual_heat_shock),
         annual_heat_shock_cumsum_lag1 = shift(annual_heat_shock_cumsum, 1),
         annual_heat_shock_cumsum_lag2 = shift(annual_heat_shock_cumsum, 2),
         seasonal_heat_shock_cumsum = cumsum(seasonal_heat_shock),
         seasonal_heat_shock_cumsum_lag1 = shift(seasonal_heat_shock_cumsum, 1),
         seasonal_heat_shock_cumsum_lag2 = shift(seasonal_heat_shock_cumsum, 2)) %>% 
  mutate(annual_heat_shock_window = RcppRoll::roll_sum(annual_heat_shock, n = 5, align = "right", fill = NA),
         annual_heat_shock_window_lag2 = shift(annual_heat_shock_window, 2),
         seasonal_heat_shock_window = RcppRoll::roll_sum(seasonal_heat_shock, n = 5, align = "right", fill = NA),
         seasonal_heat_shock_window_lag2 = shift(seasonal_heat_shock_window, 2))


#### Temperature spatial terms 

## Load in temperature series to residualize temperature in neighbours 

temp_full_series <- full_join(temp_raw, hottest_season, by = "iso3c") %>% 
  arrange(iso3c, year, month, day) %>% 
  select(iso3c, year, month, day, season, hottest_season, temperature) %>% 
  # Get annual average temperature
  group_by(iso3c, year) %>% 
  mutate(annual_average_temp = mean(temperature, na.rm=T)) %>% 
  ungroup() %>% 
  # Select only the days in the hottest season
  filter(season == hottest_season) %>% 
  # Temperature years run December--November for seasonality
  # Offset December to match seasonal year to year
  mutate(season_year = ifelse(month == 12, year + 1, year)) %>% 
  # Drop year and rename season_year for merging
  select(-year) %>% 
  rename(year = season_year) %>% 
  filter(year<2019) %>% 
  # Group by country and season_year
  group_by(iso3c, year) %>% 
  # Create the annual average seasonal temp in the hottest season
  mutate(hottest_season_temp = mean(temperature)) %>% 
  ungroup() %>% 
  # Collapse to the year level, drop days and months
  distinct(iso3c, year, .keep_all = T) %>% 
  # Arrange for lags
  arrange(iso3c, year) %>% 
  select(iso3c, year, annual_average_temp, hottest_season_temp) %>% 
  group_by(iso3c) %>% 
  # Set lags
  do(data.frame(., setNames(shift(.$annual_average_temp, 1:3),
                            paste0('annual_average_temp_lag', 1:3)))) %>% 
  do(data.frame(., setNames(shift(.$hottest_season_temp, 1:3),
                            paste0('hottest_season_temp_lag', 1:3))))


## Load in contiguity data

# Raw contiguity
contdird_raw <- read_csv("data/contdird.csv")

# Tidy
contdird_df <- contdird_raw %>% 
  filter(year == max(year, na.rm=T)) %>%
  mutate(contiguous = ifelse(conttype>0, 1, 0)) %>% 
  mutate(iso3_state1 = countrycode::countrycode(.$state1no, origin = 'cown', destination = 'iso3c'),
         iso3_state2 = countrycode::countrycode(.$state2no, origin = 'cown', destination = 'iso3c')) %>% 
  select(iso3_state1, iso3_state2, contiguous)

# Extract states list from temperature data
states_list_analysis <- temp_full_series %>% 
  filter(!is.na(annual_average_temp),
         year == 2005,
         iso3c!="TWN") %>% 
  select(iso3c) %>% 
  as.vector()

# Unique, sorted states list
states_list <- cbind.data.frame(iso3c = unique(contdird_df$iso3_state1))
states_list <- states_list %>% filter(iso3c %in% states_list_analysis$iso3c)
states_list <- sort(unique(states_list$iso3c))

# Make contiguity dyadic data
cont_df <- data.frame(iso3_state1 = rep(states_list, times = length(states_list)),
                      iso3_state2 = rep(states_list, each = length(states_list))) %>% 
  left_join(., contdird_df, by = c("iso3_state1", "iso3_state2")) %>% 
  mutate(contiguous = ifelse(is.na(contiguous), 0, contiguous)) %>% 
  arrange(iso3_state1, iso3_state2)
# cont_df %>% .[1:10,]

# Row-standardize and make a matrix
cont_matrix <- cont_df %>% 
  pivot_wider(names_from = iso3_state2, values_from = contiguous) %>% 
  column_to_rownames("iso3_state1") %>% 
  as.data.frame() %>% 
  mutate(row_total = rowSums(across(where(is.numeric)))) %>% 
  mutate_at(vars(AFG:ZWE), ~./row_total) %>% 
  select(-row_total) %>% 
  as.matrix()
diag(cont_matrix) <- 0

## Calculate each country's residual of temperature, for lags and seasons
temp_demeaned <- temp_full_series %>% 
  ungroup() %>% 
  filter(iso3c %in% states_list) %>% 
  filter(!is.na(annual_average_temp)) %>% 
  filter(year>=1990) %>% 
  # Annual average temperature
  mutate(resid_annual_lag1 = residuals(lm(annual_average_temp_lag1 ~ 
                                            as.factor(iso3c) + 
                                            as.factor(year),
                                          data = .))) %>% 
  mutate(resid_annual_lag2 = residuals(lm(annual_average_temp_lag2 ~ 
                                            as.factor(iso3c) + 
                                            as.factor(year),
                                          data = .))) %>% 
  mutate(resid_annual_lag3 = residuals(lm(annual_average_temp_lag3 ~ 
                                            as.factor(iso3c) + 
                                            as.factor(year),
                                          data = .))) %>% 
  # Hottest season average temperature
  mutate(resid_season_lag1 = residuals(lm(hottest_season_temp_lag1 ~ 
                                            as.factor(iso3c) + 
                                            as.factor(year),
                                          data = .))) %>% 
  mutate(resid_season_lag2 = residuals(lm(hottest_season_temp_lag2 ~ 
                                            as.factor(iso3c) + 
                                            as.factor(year),
                                          data = .))) %>% 
  mutate(resid_season_lag3 = residuals(lm(hottest_season_temp_lag3 ~ 
                                            as.factor(iso3c) + 
                                            as.factor(year),
                                          data = .))) %>% 
  select(iso3c, year, annual_average_temp, starts_with("resid")) %>% 
  arrange(iso3c)
# head(temp_demeaned)

### Now make the spatial weights

# Make object for loops

# To hold loop outputs
spatial_temp_annual_lag1 <-  as.data.frame(matrix(NA, nrow = 188, ncol = 28))
spatial_temp_annual_lag2 <-  as.data.frame(matrix(NA, nrow = 188, ncol = 28))
spatial_temp_annual_lag3 <-  as.data.frame(matrix(NA, nrow = 188, ncol = 28))
spatial_temp_season_lag1 <-  as.data.frame(matrix(NA, nrow = 188, ncol = 28))
spatial_temp_season_lag2 <-  as.data.frame(matrix(NA, nrow = 188, ncol = 28))
spatial_temp_season_lag3 <-  as.data.frame(matrix(NA, nrow = 188, ncol = 28))

# List of years to filter on and loop through
year_list <- unique(temp_demeaned$year)

for (yearr in year_list) {
  residuals <- temp_demeaned %>% 
    filter(year == yearr) 
  
  # Annual, lag 1
  hold_spatial_weight <- cont_matrix %*% residuals$resid_annual_lag1
  spatial_temp_annual_lag1[,yearr-1989] <- hold_spatial_weight
  # Annual, lag 2
  hold_spatial_weight <- cont_matrix %*% residuals$resid_annual_lag2
  spatial_temp_annual_lag2[,yearr-1989] <- hold_spatial_weight
  # Annual, lag 3
  hold_spatial_weight <- cont_matrix %*% residuals$resid_annual_lag3
  spatial_temp_annual_lag3[,yearr-1989] <- hold_spatial_weight
  # Season, lag 1
  hold_spatial_weight <- cont_matrix %*% residuals$resid_season_lag1
  spatial_temp_season_lag1[,yearr-1989] <- hold_spatial_weight
  # Season, lag 2
  hold_spatial_weight <- cont_matrix %*% residuals$resid_season_lag2
  spatial_temp_season_lag2[,yearr-1989] <- hold_spatial_weight
  # Annual, lag 3
  hold_spatial_weight <- cont_matrix %*% residuals$resid_season_lag3
  spatial_temp_season_lag3[,yearr-1989] <- hold_spatial_weight
}

# Merge all the spatial weights together
spatial_temp_df <- full_join(spatial_temp_annual_lag1 %>% 
                               pivot_longer(names_to = "names", values_to = "annual_W_lag1", V1:V29) %>% 
                               mutate(iso3c = rep(states_list, times = length(year_list)),
                                      year = rep(year_list, each = length(states_list))) %>% 
                               select(iso3c, year, 2),
                             spatial_temp_annual_lag2 %>% 
                               pivot_longer(names_to = "names", values_to = "annual_W_lag2", V1:V29) %>% 
                               mutate(iso3c = rep(states_list, times = length(year_list)),
                                      year = rep(year_list, each = length(states_list))) %>% 
                               select(iso3c, year, 2),
                             by = c("iso3c", "year")) %>% 
    full_join(., spatial_temp_annual_lag3 %>% 
                pivot_longer(names_to = "names", values_to = "annual_W_lag3", V1:V29) %>% 
                mutate(iso3c = rep(states_list, times = length(year_list)),
                       year = rep(year_list, each = length(states_list))) %>% 
                select(iso3c, year, 2),
              by = c("iso3c", "year")) %>% 
    full_join(., spatial_temp_season_lag1 %>% 
                pivot_longer(names_to = "names", values_to = "season_W_lag1", V1:V29) %>% 
                mutate(iso3c = rep(states_list, times = length(year_list)),
                       year = rep(year_list, each = length(states_list))) %>% 
                select(iso3c, year, 2),
            by = c("iso3c", "year")) %>% 
  full_join(., spatial_temp_season_lag2 %>% 
              pivot_longer(names_to = "names", values_to = "season_W_lag2", V1:V29) %>% 
              mutate(iso3c = rep(states_list, times = length(year_list)),
                     year = rep(year_list, each = length(states_list))) %>% 
              select(iso3c, year, 2), 
            by = c("iso3c", "year")) %>% 
  full_join(., spatial_temp_season_lag3 %>% 
              pivot_longer(names_to = "names", values_to = "season_W_lag3", V1:V29) %>% 
              mutate(iso3c = rep(states_list, times = length(year_list)),
                     year = rep(year_list, each = length(states_list))) %>% 
              select(iso3c, year, 2),
            by = c("iso3c", "year")) 

# Remove intermediate objects
rm(contdird_raw, contdird_df, temp_full_series,
   states_list_analysis, states_list, 
   cont_df, cont_matrix, temp_demeaned, 
   spatial_temp_annual_lag1, spatial_temp_annual_lag2, 
   spatial_temp_annual_lag3, spatial_temp_season_lag1, 
   spatial_temp_season_lag2, spatial_temp_season_lag3, 
   year_list, yearr, residuals, 
   hold_spatial_weight)

# Merge these spatial weights into the analysis dataset

#### Join these all together

unique_iso <- unique(gcg$iso3c)
analysis_df <- full_join(temp_df, spatial_temp_df, by = c("iso3c", "year")) %>% 
  full_join(., climate_laws, by = c("iso3c", "year")) %>% 
  full_join(., ecp, by = c("iso3c", "year")) %>%
  full_join(., wdi_renew, by = c("iso3c", "year")) %>%
  full_join(., climate_policies, by = c("iso3c", "year")) %>% 
  full_join(., cop, by = c("iso3c", "year")) %>% 
  full_join(., gcg, by = c("iso3c", "year")) %>% 
  full_join(., finance_donor, by = c("iso3c", "year")) %>% 
  full_join(., finance_recipient, by = c("iso3c", "year")) %>% 
  left_join(., democracy, by = c("iso3c")) %>% 
  left_join(., rich, by = c("iso3c")) %>% 
  filter(year>1989) %>% 
  # Subset the data to only ones with climate institution membership data 
  filter(iso3c %in% unique_iso) %>% 
  select(-carbon_priced, -scaled_carbon_priced) %>% 
  mutate_at(vars(matches("scaled")), ~as.numeric(.)) %>% 
  as.data.frame()

# Remove intermediate objects
rm(climate_laws, climate_policies,
   cop, democracy, ecp, finance_donor, finance_recipient,
   gcg, hottest_season, rich, temp_df, temp_raw,
   wdi_renew, myFiles, unique_iso)
# rm(spatial_temp_df)
# rm(list=setdiff(ls(), "analysis_df"))

